*****This estimates fully flexible diff in diffs around the cutoffss and for the whole sample

global user "diegove"
global dirdata "C:\Users\\$user\Dropbox\CCT BJP\data"
global dir1  "$dirdata\Final data"
global dir2  "$dirdata\Results\graphs"
cd "$dir2"

use "$dir1\pooled_hh.dta", clear

** Exclude the 5obs for which we observe a -5 and -6 years (8th-7th graders).
keep if  tau>-5
keep if work_m!=. & work_p!=.

g age=age_m if hhmale==0
replace age=age_p if hhmale==1
g spanish=spanish_m if hhmale==0
replace spanish=spanish_p if hhmale==1
global covs "age_m spanish_m schooling_m age_p spanish_p schooling_p hhmale nchild5 urban age_sac np c.age_sac#i.year c.np#i.year i.departamento i.departamento#i.year i.year"



egen work_heads=rowtotal(work_m work_p)
g work_one=work_heads==1
g work_two=work_heads==2

areg  work_one tau_1 tau_2  tau_3 tau_4 tau_5 tau_7 tau_8 tau_9 tau_10 i.min_schooling   $covs  [w=factor]  if (min_schooling>=1 & min_schooling<9) , a(departamento) cluster(municipio)
mat V=e(V)
preserve
clear
set obs 10
gen coef=.
replace coef=_b[tau_1] in 1
replace coef=_b[tau_2] in 2
replace coef=_b[tau_3] in 3
replace coef=_b[tau_4] in 4
replace coef=_b[tau_5] in 5
replace coef=0 in 6
replace coef=_b[tau_7] in 7
replace coef=_b[tau_8] in 8
replace coef=_b[tau_9] in 9
replace coef=_b[tau_10] in 10
gen sd2=.

replace sd2=V[1,1] in 1
replace sd2=V[2,2] in 2
replace sd2=V[3,3] in 3
replace sd2=V[4,4] in 4
replace sd2=V[5,5] in 5

replace sd2=0 in 6
replace sd2=V[6,6] in 7
replace sd2=V[7,7] in 8
replace sd2=V[8,8] in 9
replace sd2=V[9,9] in 10

gen year=.
replace year=-6 in 1
replace year=-5 in 2
replace year=-4 in 3
replace year=-3 in 4
replace year=-2 in 5
replace year=-1 in 6
replace year=0 in 7 
replace year=1 in 8 
replace year=2 in 9
replace year=3 in 10
g ub=coef+1.96*sqrt(sd2)
g lb=coef-1.96*sqrt(sd2)
tw scatter coef year  if year>-5, sort || rcap ub lb year  if year>-5, sort  yline(0) ytitle(Prob.) legend(off)  xtitle("")   title("A-Only one parent works")   ylabel(#10) ylabel(,labsize(small)) xlabel(#9) graphregion(color(white)) xline(-1,lcolor(maroon) lpattern(dash))
graph export Figure6_top.eps , replace
restore


areg  work_two tau_1 tau_2  tau_3 tau_4 tau_5 tau_7 tau_8 tau_9 tau_10 i.min_schooling   $covs  [w=factor]  if (min_schooling>=1 & min_schooling<9) , a(departamento) cluster(municipio)
mat V=e(V)
preserve
clear
set obs 10
gen coef=.
replace coef=_b[tau_1] in 1
replace coef=_b[tau_2] in 2
replace coef=_b[tau_3] in 3
replace coef=_b[tau_4] in 4
replace coef=_b[tau_5] in 5
replace coef=0 in 6
replace coef=_b[tau_7] in 7
replace coef=_b[tau_8] in 8
replace coef=_b[tau_9] in 9
replace coef=_b[tau_10] in 10
gen sd2=.

replace sd2=V[1,1] in 1
replace sd2=V[2,2] in 2
replace sd2=V[3,3] in 3
replace sd2=V[4,4] in 4
replace sd2=V[5,5] in 5

replace sd2=0 in 6
replace sd2=V[6,6] in 7
replace sd2=V[7,7] in 8
replace sd2=V[8,8] in 9
replace sd2=V[9,9] in 10

gen year=.
replace year=-6 in 1
replace year=-5 in 2
replace year=-4 in 3
replace year=-3 in 4
replace year=-2 in 5
replace year=-1 in 6
replace year=0 in 7 
replace year=1 in 8 
replace year=2 in 9
replace year=3 in 10
g ub=coef+1.96*sqrt(sd2)
g lb=coef-1.96*sqrt(sd2)
tw scatter coef year  if year>-5, sort || rcap ub lb year  if year>-5, sort  yline(0) ytitle(Prob.) legend(order(1 "Coefficient" 2 "95 % CI"))  xtitle("Time to treatment")   title("B-Both parents work")  ylabel(#10) ylabel(,labsize(small)) xlabel(#9) graphregion(color(white)) xline(-1,lcolor(maroon) lpattern(dash))
graph export Figure6_bottom.eps , replace
restore
